subroutine Wsigma_GS(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,&
    di,f,l1,v1,v2,v3,v4,v5,ggl2,ggu2,di2,vertex3D,nVer3D,nInCat3D,igT3D,jgT3D,&
    ggT3D,log_txt,numArea,zLayers,nLayers,vnRZ,wnRZ,wnXYZ)
    common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,m2D,mTerm2d,nodesNotTermCount2D
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,edgesNotTermCount3D,&
    edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    integer elemsCount2D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D
    integer typeTask,calcWK,calcVP,vn2d_dat,vn3d_dat,log_txt
    real*8 vertex2D,vertex3D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,wnRZ,wnXYZ
    real*8 eps2d,eps3d,eps3dNS,radiusRED,ggT3D,vnRZ,zLayers
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(2,nMaterials)
    dimension igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),f(nodesNotTermCount2D),l1(nodesBoundCount2D),&
    v1(nodesNotTermCount2D),v2(nodesNotTermCount2D),v3(nodesNotTermCount2D),v4(nodesNotTermCount2D),&
    v5(nodesNotTermCount2D),ggl2(m2d),ggu2(m2d),di2(nodesNotTermCount2D),wnRZ(nodesCount2D),&
    vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),nInCat3D(elemsCount3D),wnXYZ(nodesCount3D),&
    igT3D(nodesTermCount3D+1),jgT3D(mTerm3d),ggT3D(mTerm3d),vnRZ(nodesCount2D),zLayers(nLayers)

    call clearVectorDouble(wnRZ,nodesCount2D)
    call clearVectorDouble(wnXYZ,nodesCount3D)

    !Assembling global matrix
    call globalForWsigma_GS(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,&
    params,ggT,igT,jgT,vnRZ,numArea,zLayers,nLayers)

    !Taking into account boundary condition
    call firstCondition(di,f,l1,nodesNotTermCount2D,nodesBoundCount2D,1)
    
    !Solving SLAE
    call LOS3(wnRZ,f,ig,jg,ggl,ggu,di,nodesNotTermCount2D,m2d,&
    maxiter,eps2d,v1,v2,v3,v4,ggl2,ggu2,di2,v5,0,iter,residual)
    call calculateValueForNotTermNodes(wnRZ,igT,jgT,ggT,nodesCount2D,&
    mTerm2d,nodesTermCount2D,nodesNotTermCount2D)

    !Recalculating solution to XYZ
    call v3d(vertex3D,nVer3D,vertex2D,nVer2D,wnRZ,wnXYZ)
    call calculateValueForNotTermNodes(wnXYZ,igT3D,jgT3D,ggT3D,nodesCount3D,&
    mTerm3d,nodesTermCount3D,nodesNotTermCount3D)

2   format(A,i5,' residual =',E)
    write(*,2) 'VP: end W_n:       iter =',iter,residual
    write(log_txt,2)'VP: end W_n:       iter =',iter,residual
end
!*************************************************************************
subroutine globalForWsigma_GS(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,&
    params,ggT,igT,jgT,vnRZ,numArea,zLayers,nLayers)
    common /dim2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,m2D,mTerm2d,nodesNotTermCount2D
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,edgesNotTermCount3D,&
    edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    integer elemsCount2D,edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,edgeTerm2d
    integer elem1,elem2,calcWK,calcVP,typeTask
    real*8 vertex2D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,wnRZ
    real*8 g,r1,r2,z1,z2,sigma,b,eps2d,eps3d,eps3dNS,radiusRED,v,vnRZ,zLayers
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(2,nMaterials),v(4)
    dimension g(4),b(4,4),igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),&
    ig(nodesNotTermCount2D+1),jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),&
    f(nodesNotTermCount2D),vnRZ(nodesCount2D),zLayers(nLayers)

    call clear(ggl,ggu,di,f,nodesNotTermCount2D,m2d)
    do i=1,elemsCount2D
        r1=vertex2D(1,nVer2D(1,i))
        r2=vertex2D(1,nVer2D(2,i))
        z1=vertex2D(2,nVer2D(1,i))
        z2=vertex2D(2,nVer2D(3,i))
        sigma=params(2,nInCat2D(i))

        if(r1.eq.0.d0)r1=1d-3

        do j=1,4
            v(j)=vnRZ(nVer2D(j,i))
        enddo

        !Assembling local matrix
        call LocalMatrixForWsigma_GS(b,g,r1,r2,z1,z2,sigma,v,numArea,zLayers,nLayers)
        !Adding local matrix to global
        call AddLocal(di,ig,jg,ggl,ggu,f,b,g,nVer2D(1,i),nodesNotTermCount2D,&
        m2d,4,igT,jgT,ggT,nodesTermCount2D,mTerm2d,1)
    enddo
end
!*************************************************************************
subroutine LocalMatrixForWsigma_GS(b,g,r1,r2,z1,z2,sigma,v,numArea,&
    zLayers,nLayers)
    dimension b(4,4),br(2,2),cr(2,2),bz(2,2),cz(2,2),v(4),g(4),zLayers(nLayers)
    real*8 b,r1,r2,z1,z2,br,cr,bz,cz,sigma,v,g,zLayers
    integer delta,d

    call local_bc_r(br,cr,r1,r2)
    call local_bc_z(bz,cz,z1,z2)

    do i=1,4
        do j=1,4
          mi=mu(i)
          mj=mu(j)
          ni=nu(i)
          nj=nu(j)
          b(i,j)=br(mi,mj)*cz(ni,nj)+bz(ni,nj)*cr(mi,mj)
        enddo
    enddo

    call mult_sq(g,b,v,4)
    d = delta(numArea,z1,z2,zLayers,nLayers)
    do i=1,4
        g(i)=g(i)*sigma*d
        do j=1,4
            b(i,j)=b(i,j)*sigma
        enddo
    enddo
end
!*************************************************************************
!*************************************************************************
!*************************************************************************